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Abstract 

The sagebrush lizards ( Sceloporus graciosus group) consist of four taxa ( S. graciosus graciosus, S. graciosus gracilis, S. 
graciosus vandenburgianus , and S. arenicolus ) distributed in western North America. Of these, S. arenicolus is morpho¬ 
logically, behaviorally, and ecologically distinct as well as geographically disjunct from the other taxa, occurring only in 
the Mescalero-Monahans Sandhills of southeastern New Mexico and adjacent Texas. Sceloporus arenicolus is a taxon of 
concern because of its small range and habitat alteration due to land use practices. Understanding evolutionary relation¬ 
ships among members of the S. graciosus group, and especially S. arenicolus, has important implications for conservation. 
We examine the phylogenetic relationship of S. arenicolus relative to the three recognized subspecies of S. graciosus at 
mitochondrial and nuclear loci for populations sampled throughout the ranges of these taxa. Additionally, we estimate the 
divergence time and clade age of S. arenicolus. We find that the S. graciosus group is in need of major taxonomic revision, 
and also confirm that S. arenicolus is a genetically distinct and divergent lineage. These results bear important conse¬ 
quences for conservation and management. 

Key words: Sceloporus graciosus, Mescalero-Monahans Sand Dunes, Phrynosomatidae, shinnery oak, evolutionarily 
significant unit (ESU), sagebrush lizard 


Introduction 

The sagebrush lizards are comprised of two species in the graciosus group of the phrynosomatid genus Sceloporus 
Wiegmann (1828)— Sceloporus graciosus Baird and Girard (1852) containing three subspecies S. g. graciosus 
Baird and Girard (1852), S. g. gracilis Baird and Girard (1852), and S. g. vandenburgianus Cope (1896) and S. 
arenicolus Degenhardt and Jones (1972). The taxonomy of this group has fluctuated considerably; Sceloporus g. 
vandenburgianus has previously been considered a distinct species from S. graciosus (Collins 1991), and likewise, 
S. arenicolus has previously been considered a geographical variant and subspecies of S. graciosus (Degenhardt & 
Jones 1972). As currently defined, Sceloporus graciosus occurs through most of the Great Basin of the western 
USA, to the coast of California, to northern Arizona and northwestern New Mexico (Figure 1; Stebbins 2003; Ryan 
2009). The geographic range of Sceloporus arenicolus is disjunct from the ranges of the S. graciosus subspecies 
and occurs in southeastern New Mexico and adjacent Texas (Stebbins 2003; Fitzgerald & Painter 2009; Laurencio 
& Fitzgerald 2010). 

The morphological and ecological distinctiveness of S. arenicolus from S. graciosus is well-studied. 
Specimens of S. arenicolus were originally identified by Sabath (1960), who reported it from Texas and New 
Mexico as a range extension of S. graciosus with distinct coloration patterns and femoral pore counts. Kerfoot 
(1968) conducted a comprehensive study of morphological variation among 1,311 specimens from throughout the 
eastern portion of the range of S. graciosus including the “southeastern sand dune isolates” which were populations 
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of S. arenicolus, then considered S. g. graciosus. Kerfoot’s (1968) analyses identified S. arenicolus as a clear 
outlier from all other populations of S. graciosus based on distinct coloration and scalation of the body and head. 
Kerfoot (1968) noted, “Dorsal scale and midbody scale counts of these lizards [5. arenicolus ] are unusually low, 
whereas the total number of femoral pores is unusually high.” These populations were formally described as a 
subspecies, S. g. arenicolus, by Degenhardt and Jones (1972) based on morphology and coloration and Cole (1975) 
used karyotype data to confirm that S. arenicolus was more closely related to S. graciosus than S. undulatus (now 
S. consobrinus), a sympatric species. Based on the criterion of allopatry, Collins (1991) elevated three subspecies 
of the S. graciosus group ( S. g. arenicolus, S. g. graciosus, S. g. vandenburgianus) to the level of species. Although 
this criterion was arbitrarily applied and a formal analysis was lacking, subsequent works have treated S. arenicolus 
as a valid species (Smith et al. 1992; Degenhardt et al. 1996; Fitzgerald & Painter 2009; Ryan 2009; deQuieroz & 
Reeder 2012). 



FIGURE 1. Collection localities for samples from the Sceloporus graciosus group included in this study. Colored symbols 
correspond to clade membership on Figure 3; three individuals for which we only have sequence data atR55 are indicated with 
stars. Putative species and subspecies boundaries are shaded for the members of the Sceloporus graciosus group. Sceloporus 
graciosus graciosus'. dark gray distribution, S. g. gracilis', light gray distribution, S. g. vandenburgianus'. brown distribution, 
and S. arenicolus : black distribution. 
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Despite the apparent ecological and morphological distinctiveness of S. arenicolus and presumably long 
separation from other members of the S. graciosus group, the phylogenetic placement of 5. arenicolus has mainly 
been explored within the context of broader systematic studies of the genus (Wiens & Reeder 1997; Leache 2010; 
Wiens et al. 2010). A single study by Frabotta (2002) examined the mitochondrial relationships among populations 
from all members of the S. graciosus group and suggested that S. arenicolus evolved from within S. graciosus, 
though sampling from eastern S. g. graciosus and S. arenicolus was limited. 

Further investigation of these proposed relationships with greater population level sampling throughout the 
geographic ranges of each group is warranted both for basic and applied reasons. Better understanding of the 
phylogenetic placement of S. arenicolus would clarify its relationships to other members of the S. graciosus group 
and give finer resolution to the clade within the phylogeny of Sceloporus. Additionally, clear understanding of the 
phylogenetic and taxonomic position of S. arenicolus bears important consequences for conservation. Sceloporus 
arenicolus is endemic to the Mescalero-Monahans Sandhills, a small and restricted ecosystem. Within this 
ecosystem S. arenicolus is an ecological specialist that only occurs in dunes stabilized by growing shinnery oak 
('Quercus havardii) that contain open sandy depressions (Fitzgerald & Painter 2009; Ryberg et al. 2012). Land use 
in the region has contributed to loss and fragmentation of shinnery dunes habitat, and the conservation status of the 
species is of growing concern to natural resource agencies, oil and gas interests, and landowners (Smolensky & 
Fitzgerald, 2010; 2011). Sceloporus arenicolus has been of concern to the US Fish and Wildlife Service since 1982, 
and it was proposed for listing as Endangered in December 2010 (U.S. Fish and Wildlife 2010). An understanding 
of the uniqueness, history, and divergence of S. arenicolus based on phylogenetic analyses will help inform and 
justify efforts to conserve it. 

The phylogenetic relationships at the scale of species groups within the genus Sceloporus are well understood 
(Wiens & Reeder 1997; Leache 2010; Wiens et al. 2010) and information also exists on population differentiation 
within S. arenicolus (Chan et al. 2009). As such, it is not our aim to re-evaluate deeper relationships within 
Sceloporus, nor fmer-scale patterns of phylogenetic relationships within S. arenicolus. The purpose of this study is 
to use multilocus sequence data to clarify the phylogenetic relationship of S. arenicolus relative to populations of S. 
graciosus. In addition, we estimate the clade age of extant S. arenicolus and the divergence time from its most 
recent common ancestor (MRCA) to evaluate the evolutionary history of this group. This work will enhance our 
understanding of relationships in the S. graciosus group as well as provide the first comprehensive phylogenetic 
analysis of the group that clarifies the position of S. arenicolus. 


Methods 

Sampling and molecular genetic data collection. Tissues samples for S. graciosus and S. arenicolus were 
obtained from vouchered specimens in museum and personal collections and from individuals captured in the field 
(toe and tail preserved in 95% EtOH) (Appendix 1). Whole genomic DNA was extracted from liver samples and 
tail and toe clips using the DNeasy Blood and Tissue Kit (Qiagen). Five protein-coding gene regions were targeted 
for sequencing using PCR amplification. These included two mitochondrial regions, NADH dehydrogenase 1 
( ND1 ) with the primers tMet and 16dR (Leache & Reeder 2002) and cytochrome-b ( cyt-b ) using the primers 
L14724 and HI5915 (Irwin et al. 1991), in addition to three nuclear regions, the pinin gene ( PNN ), recombination 
activating gene-1 ( RAG1 ), and RNA fingerprint protein 35 ( R35 ) (primers from Leache 2010). 

PCR was done in 15 pi reactions with 1 x Buffer, 2.0 mM MgCl„ 0.33 mM each dNTP (0.53 mM for RAG1), 
0.4-0.6 pM each primer, and 0.375 U (mtDNA), 0.5 U ( R35 and PNN), or 0.75 U ( RAG1 ) Taq polymerase (Sigma 
JumpStart). Reactions for RAG1 and PNN additionally included 1% DMSO. Amplification cycles consisted of an 
initial denaturation step at 94°C for 3 minutes followed by 35 cycles (40 for RAG1) of 94°C for 1 min (30 sec for 
nuclear genes), locus specific annealing temperature (54°C ND1, 45°C cyt-b, 52°C RAG1, 58°C R35 and PNN) for 
1 min (30 sec for nuclear genes) and a 72°C extension for 90 sec. Cycles were followed by a final extension of 
72°C for 10 min. 

PCR product was cleaned with 0.5 pL of ExoSapIT (USB) and sequenced in 1/16 th reactions using the same 
primers used for amplification. Sequence data were collected on an ABI 3730x1. Chromatograms were assembled 
and checked for errors in Geneious 5.6 (Biomatters). We included sequence data from both mitochondrial genes for 
three phrynosomatid relatives ( Urosaurus, Uta, and Phrynosoma) and three Sceloporus species {S. jarrovii, S. 
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merriami, and S. occidentalis) representing major clades in the analyses ofLeache (2010) (Appendix 2). We used 
the MAFFT plug-in in Geneious to align sequences for each gene region. Alignments were adjusted by hand and 
checked to ensure the absence of premature stop-codons; new sequences were deposited in GenBank (KC853767 - 
KC854153). For each nuclear locus, we constructed parsimony haplotype networks in TCS 1.2 (Clement et al. 
2000) using a connection limit of 95%. 

Phylogenetic reconstruction. Phylogenetic analyses were conducted on the concatenated mitochondrial DNA 
dataset under both Bayesian and Maximum Likelihood (ML) statistical frameworks. Redundant mitochondrial 
haplotypes were omitted from the aligmnent prior to analysis to reduce computational time. We used DT-ModSel 
(Minin et al. 2003) to detennine the best-fit model of DNA sequence evolution for each codon position of each 
gene. Partitioned Bayesian phylogenetic analyses were conducted in MrBayes 3.2.1 (Ronquist & Fluelsenbeck 
2003). Two independent final runs each consisted of a single chain run sampled every 1,000 generations for 25 
million generations. We confirmed adequate mixing within runs and convergence between runs in Tracer 1.5 
(Rambaut & Drummond 2007) and summarized the final 75% posterior distribution of trees on the half-compatible 
consensus phylogeny. One-hundred partitioned ML bootstrap replicates were conducted in Garli 2.0 (Zwickl 
2006). Each bootstrap replicate consisted of five independent search replicates with search tennination thresholds 
of 5,000 generations and 0.001 LnL units. Bootstrap scores were summarized on the consensus phylogeny from the 
Bayesian phylogenetic analysis using the SumTree script from the DendroPy package (Sukumaran & Flolder 
2010 ). 

Divergence time estimation. We additionally estimated the mitochondrial divergence times among all 
individuals sampled from the S. graciosus group in BEAST 1.7.5 (Drummond & Rambaut 2007). We assumed 
uncorrelated relaxed lognormal clocks for separate ND1 and cyt-b partitions each with normally distributed prior 
substitution rates (mean of 1.5% per myr, standard deviation of 0.5% myr). This brackets mitochondrial 
substitution rates reported for other iguanian lizards (Bryson et al. 2011; Chan et al. 2012). We conducted analyses 
under a Yule speciation tree model. Replicate analyses were conducted to ensure adequate mixing and convergence 
and the final run consisted of 50 million generations sampled every 5,000 generations. We summarized the clade 
ages after discarding the first 10% of samples as bumin. 


Results 

We recovered the entire coding region of ND1 (990 base pairs, bp) for 85 individuals from the S. graciosus group 
and partial cyt-b sequences (1128 bp) for 81 individuals. Of the 2,118 bp of mtDNA, 748 sites were variable and 
617 were parsimony infonnative. Nuclear regions were much less variable: for 930 bp of RAW for 64 individuals, 
nine sites were variable and four were parsimony informative. Eighty-eight individuals sequenced for R35 had 658 
bp; of these 25 were variable and 11 were parsimony informative. Finally, we recovered 1,046 bp for RAG1 from 
70 individuals of which 14 of 29 variable sites were parsimony informative. 

Flaplotype networks for the nuclear loci showed that common haplotypes were shared among members of the 
S. graciosus group (Figure 2). Sceloporus arenicolus had unique alleles at all three loci, but shared a common allele 
with all subspecies of S. graciosus at PNN and with S. g. graciosus at RAG1. 

We partitioned the mtDNA sequence data by gene and codon position for phylogenetic analyses. For cyt b, an 
FIKY+I+T model was applied to first and second codon positions, and a TrN+ T model was applied to the third 
position. For ND1, TVM+ T, F1KY+I, and TIM+I+ T were applied to first, second, and third codon positions, 
respectively. Bayesian and ML analyses were highly concordant (Figure 3). We found strong support for the 
monophyly of the S. graciosus group as well as the monophyly of S. arenicolus (PPM, Bootstrap = 100). None of 
the other putative taxa ( S. g. graciosus, S. g. gracilis, and S. g. vandenburgianus ) were supported as monophyletic. 

Divergence time analyses assuming a mean substitution rate of 1.5% per mya estimated the mean root age of 
the S. graciosus group to be 9.21 million years (Ma) with a 95% highest posterior density (HPD) of 5.70 - 13.91 
Ma. The mean divergence time between S. arenicolus and the MRCA was 2.55 Ma (95% HPD: 1.43 - 3.93 Ma) 
while the mean clade age of S. arenicolus was 0.66 Ma (95% HPD: 0.34 - 1.06 Ma). 
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FIGURE 2. Minimum spanning haplotype networks for all S. graciosus group samples sequenced at each of three nuclear loci. 
Size of each circle corresponds to the frequency of that haplotype. Shading corresponds to clade membership in Figure 3. 


Discussion 

Consistent with previous studies focusing on broader systematic relationships (Wiens & Reeder 1997; Leache 
2010; Wiens et al. 2010), we find strong support for a monophyletic S. graciosus group. We additionally find 
unequivocal support for the monophyly of S. arenicolus at mtDNA loci. In contrast, S. g. graciosus, S. g. gracilis, 
and S. g. vandenburgianus are each recovered as paraphyletic or polyphyletic indicating that the S. graciosus group 
as a whole requires major taxonomic revision. The only other phylogenetic study of the S. graciosus group to 
include population sampling to date (Frabotta 2002) suggested that S. arenicolus fell within S. g. graciosus. That 
study had limited samples of eastern S. g. graciosus and in this study we were able to incorporate numerous 
samples from the eastern range of S. g. graciosus including samples from the populations in closest proximity to S. 
arenicolus. As such, our results corroborate and enhance previous work with greater genetic and geographic 
sampling. In both Bayesian and ML phylogenetic reconstruction of mitochondrial sequence data, S. arenicolus falls 
within a well-supported clade of eastern S. g. graciosus (Figure 3). Low nucleotide diversity at nuclear loci 
prevents phylogenetic reconstruction, but haplotype networks show that S. arenicolus does share haplotypes with 
other described taxa at two of the three nuclear gene regions examined (Figure 2). In both instances, S. arenicolus 
shares a single common haplotype at that locus and all other nuclear alleles found among S. arenicolus are unique 
to that group. Shared alleles at these two protein-coding nuclear loci are consistent with the general expectations of 
slowly evolving nuclear genes given recent divergence of S. arenicolus from within the S. graciosus group, in 
particular populations in the southeastern portion of the distribution of S. g. graciosus. 

Our phylogenetic results corroborate discussion of S. gracious group systematics by Wiens and Reeder (1997) 
stating that “Sceloporus arenicolus is clearly allopatric and diagnosable, whereas S. vandenburgianus may be 
distinct but is morphologically similar and geographically close to S. graciosus (Censky 1986).” The 
morphological distinctiveness of S. arenicolus is well-known (Sabath 1960; Kerfoot 1968; Degenhardt & Jones 
1972) and it is specialized, both ecologically and behaviorally, to live in the shinnery oak dune system where it is 
endemic. Based on a mean substitution rate of 1.5% per mya and a relaxed molecular clock, we estimate that 
divergence of S. arenicolus from its MRCA occurred approximately 2.55 Ma and extant populations coalesce to 
about 0.66 Ma. Thus, S. arenicolus has long been isolated from other members of the S. graciosus group. Given 
that it is geographically disjunct and locally endemic, it is likely this taxon is an example of divergence through 
peripheral isolation. 

Our results show S. arenicolus is a unique evolutionary lineage that is genetically distinct, in addition to being 
ecologically, morphologically, and geographically distinct. We find complex relationships among the currently 
delineated subspecies of S. graciosus with clades largely consisting of geographic clusters with no support for the 
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monophyly of any subspecies. This indicates the need for taxonomic reassessment and revision of this group and 
while such revisions are beyond the scope of this current study, we suggest that any taxonomic revisions take into 
account the species-level distinctiveness of S. arenicolus, despite its rendering S. graciosus paraphyletic (Kizirian 
& Donnelly 2004; Rieppel 2010). 
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FIGURE 3. Consensus tree from Bayesian phylogenetic analysis. Bayesian posterior probabilities (PP) and ML bootstrap 
support (BS) are noted at nodes with high support (> 95% PP and 75 BS) and at basal nodes. Nodal support of PP = 1 and BS = 
100 are indicated with a solid circle at the branch. Shaded symbols correspond to collection localities on Figure 1. 
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Conservation of S. arenicolus begins with proper identification of S. arenicolus and a stable taxonomy 
(Morrison et al. 2009). The recent proposal by the U. S. Fish and Wildlife Service to list S. arenicolus as federally 
endangered led to public commentary and confusion among stakeholder groups about the accurate taxonomic 
identity of S. arenicolus as a species (U.S. Fish & Wildlife Service 2010; 2012). As such, this study was partly 
motivated to provide a clear and stabilizing answer that would inform policy and conservation interventions (May 
1990; Funk et al. 2002). The phylogenetic evidence presented here demonstrates S. arenicolus constitutes a long- 
isolated, monophyletic lineage easily distinguished from all other populations of S. graciosus. 
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